A new approach to partial synchronization in globally coupled rotators 



O 

o 

(N 



C/3 



s 

I 



X 



p. K. MohantjE and Antonio PolitH 
Max-Planck Institute fur Physik Komplexer Systeme, 
Nothnitzer Str. 38, D-01187 Dresden, Germany 
(Dated: February 2, 2008) 

We develop a formalism to analyze the behaviour of pulse-coupled identical phase oscillators with 
a specific attention devoted to the onset of partial synchronization. The method, which allows 
describing the dynamics both at the microscopic and macroscopic level, is introduced in a general 
context, but then the application to the dynamics of leaky integrate-and-fire (LIF) neurons is anal- 
ysed. As a result, we derive a set of delayed equations describing exactly the LIF behaviour in the 
thermodynamic limit. We also investigate the weak coupling regime by means of a perturbative 
analysis, which reveals that the evolution rule reduces to a set of ordinary differential equations. 
Robustness and generality of the partial synchronization regime is finally tested both by adding 
noise and considering difi'erent force fields. 



I. INTRODUCTION 



Understanding the behavior of networks of dynamical units is a very important subject of research in various 
contexts such as information processing in the brain, or metabolic systems. Given the richness of the observed 
phenomenology and the relevance of many different ingredients (e.g., the topology of the connections, the presence 
of disorder and of delayed interactions), it is very instructive to consider even simple setups in the perspective of 
_ identifying the hopefully few mechanisms responsible for the main general phenomena. Here, we focus our attention 
^ ' on a type of collective behavior arising in globally coupled identical systems. The first evidence of collective behavior 
O ' dates back to the 70's, when it was proven that a nonzero mean field may spontaneously arise in an ensemble of 
, particles stochastically moving in a bistable potential P]. Later, it was numerically shown that macroscopic periodic 

dynamics appears in coupled stochastic Q as well as chaotic (Rossler) oscillators 3], while the first experimental 
Cn| . evidence was found in Josephson-junction arrays 0. Collective behaviour may also arise in the presence of quenched 
K*" ' disorder as revealed by the Kuramoto model [3], introduced to explain the onset of synchronization in an ensemble 
of different units Over the years it has been discovered that collective dynamics can be also chaotic, as found in 
globally coupled maps J} and oscillators ^ . 
' Here we investigate "partial synchronization" (pS), a little known phenomenon discovered in pulse-coupled, leaky 
integrate-and-fire (LIF) neurons W\. In this regime, the mean field exhibits a periodic dynamics which, at variance with 
other contexts, arises in the absence of any synchronization between the single units which behave quasi-periodically. 
, pS arises from the destabilization of a regime characterized by a constant mean field and a periodic behaviour of 
■ the single elements, their phases being equispaced. This "homogeneous regime" has been found in many different 
contexts such as Josephson devices 10], multi-mode lasers 11] and electronic circuits p^ . By following Ref. [T^ . it 
^ is now called "splay state" and its microscopic stability properties have been studied quite in detail in p^ . Since 
\. _ splay states are rather general, it is reasonable to expect that the transition to pS can be observed in a wide class of 
systems too. In fact, the preliminary evidence of pS found in Hindmarsh-Rose neurons and van der Pol oscillators 
confirms such expectations and poses the question of understanding the general conditions for the onset of pS. 
^ . In order to unravel this point, we develop a new formalism that allows treating pulse coupled oscillators and 
ILj^ ' establishes a possible basis for an extension to different coupling schemes. The approach is introduced in a general 
framework to show the potentiality of the method, but the resulting delayed equations are analyzed only in the 
specific case of LIF neurons, since in this case a simple and explicit expression is available for one of the relevant 
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?H \ variables: the "finite-time" Lyapunov exponent A. Furthermore, in the attempt of identifying the minimal model 
. . 1 able to describe the onset of pS, we have explored both analytically and numerically the weak coupling limit, showing 
that a perturbative expansion allows reducing the delayed to ordinary equations. The persistence of the transition 
to pS for an arbitrarily small coupling strength paves the way towards more general setups. A further promising 
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indication is given by the numerical observation that pS arises also for continuous force fields, which again suggests 
that the phenomenon is more general than initially conjectured @. Finally, we have added noise to the single neuron 
dynamics, finding that the macroscopic oscillations persist, though with a smaller amplitude. 



II. MODELLING IDENTICAL ROTATORS 



The investigation of globally coupled systems requires setting up a suitable (nonlinear) self-consistent approach. A 

much used approach to describe the collective behaviour of neural networks is dynamical mean field El 111 El m. 

Here, in the absence of noise and of microscopic chaos, there is no justification a priori for coarse-graining and indeed 
the transition to pS can be observed already in ensembles of finitely many oscillators ^ and can thus also be viewed 
as a bifurcation. It is therefore desirable to develop a formalism to capture both micro- and macro-scopic features. We 
start from the crossing-times as they are the most suitable variables to characterize splay states. Our approach can, in 
fact, be viewed as an extension of the method followed in Ref.|2l|, where pseudo-crossing times have been introduced 
with reference to the steady dynamics. On the other hand, it is fair to recognize that the final model equations are 
equivalent to those obtained by implementing a dynamical mean field approach, when the latter is applicable. This 
implies that in the present context one can smoothly go from the microscopic to the macroscopic level. 

We consider an ensemble of N phase oscillators (i.e. rotators), each one characterized by the phase (pi, i = 0, . . . N—1 
which evolves according to the equation 

<P, = F{^,,E{t)) (1) 

where E is a suitable mean field whose dynamics depends homogeneously on all the (pi in a way that is not crucial 
for the present discussion. The force field F{(p,E) is assumed to be positive and periodic in (p with period 1, but 
it can be discontinuous as, e.g., a sawtooth function (this is indeed the case considered in the next section). We 
start defining a suitable Poincare section, by introducing a threshold (p, which, without loss of generality, is set in 
(p = 1. Let then tm and i(m) denote, respectively, the mth crossing time (m = 1, . . . ,oo) of the threshold and the 
corresponding rotator crossin g th e threshold. Once the labels i are ordered in such a way that 0^+1 < pi, it follows 
that i{m + N) = i(m)mod N [2J| i.e., the same oscillator crosses the threshold at time tm and tm+N- 

The core of the approach consists in deriving an equation linking 5tm = im+i — tm with Stm-N in the limit of large 
N. Simple arguments show that, up to higher-order corrections in 1/iV, 

. _ 5(pi(m){tvn) ^^^{rn)(.tnl-N) _^j^i^^E})Trr, _ X-^ f (0, E{tm-N)) „A({g})r^ /r,^ 

F{l,E{tm)) F{l,E{tm)) " F{l,E{tm)) ^ ' 

where 6(j)j = <pj — (pj+i (here j is meant modulo N), while Tm — tm — tm~N is the interval between two consecutive 
threshold crossings of the same oscillator; finally, A{{E}) is the finite-time Lyapunov exponent for a trajectory starting 
at time tm~N in (p — and ending at time t„j in = 1, and thus, in general, depends on the dynamics of the field 
E. The first and last equalities in the above equation state simply that the temporal separation is equal to the phase 
separation divided by the instantaneous velocity. The middle equality follows from the observation that, being S(pi(^m) 
a small quantity, it evolves according to the linearized dynamics. In view of the possible existence of a discontinuity 
in the force field, we carefully distinguish between the (p value just before (0=1) and after {cp = 0) the threshold 
crossing. 

If the crossing-time distribution is smooth (this assumption can be verified a posteriori by integrating the resulting 
model), then 5tm has to be of the order of 1/iV and one can accordingly introduce the "instantaneous" flux 

At) ^ (3) 

Otm 

where we have dropped the subscript m in the l.h.s. as in the limit N — > oo, the time variable becomes again 
continuous. Accordingly, Eq. (jSJ simplifies to 

nt) ^)F{0,E{t-T)f ' ^' 

where the time interval T is self-consistently determined from the integral expression 

t 

dt'irit') = 1 , (5) 

t-T(t) 
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expressing the condition that from time t — T{t) to time t, all oscillators cross the threshold. An equivalent and more 
appealing equation is obtained by evaluating the time derivative of Eq. Q , namely 

In order to complete the derivation of a closed set of equations, it is necessary to include the mean field dynamics. 
In the context of pulse coupled systems, this can be easily done, since E is, by definition, the linear superposition of 
the pulses emitted by the single oscillators whenever they reach a threshold that, without loss of generality, can be 
assumed to coincide with 0. Accordingly, the mean field satisfies a linear differential equation of the type 

C(E,E,E,...)=n{t), (7) 

where the firing rate 7r(i) is nothing but the previously introduced instantaneous flux, while the structure of the 
operator C can be determined by imposing that the corresponding Green function coincides with the assumed shape 
of the emitted pulses. This completes the derivation of the final equations that are given by Eqs. H4I6I7|) . The 
model combines properties of discrete (Eq. and continuous (Eqs. HtjlTfl ) time systems, besides involving a time- 
dependent self-adjusting delay. Except for the first property, its structure resembles that of threshold delay equations 
|2l| introduced, e.g., in epidemiology and immunology (see |22j| for a similar model). 



III. AN EXAMPLE: LIF NEURONS 

In this section we specifically refer to the ensemble of LIF neurons studied already in Ref. 0. In this case, 4> 
corresponds to the (adimensional) membrane potential which evolves according to the force field 

F{4),E)=a + \{l-4))+gE{t), (8) 

where g is the coupling constant, while a and A determine the dependence of the force field on the potential (j). If (p 
reaches the threshold cf) — 1 aX time a- pulse p{t — to) is emitted (and received by all neurons), while the potential is 
reset to 0. The resetting makes it legitimate to interpret </> as a phase variable, since one can formally extend (p to the 
whole real axis by identifying (j) with + 1 and thereby see the evolution as a continuous process. Paying attention 
to the discontinuities present in the force field, Eq. writes as 

""^'^^^^'-^K + X + gEit-Tf ' 

where we have also taken into account that, because of the linearity of the force field, the Lyapunov exponent A 
appearing in Eq. is equal to —A independently of the E dynamics. 

As for the mean field, if we assume that the single pulse shape is p{t — to) = a^{t — to)e~"'-*~*°-', one finds that E 
satisfies the equation Q 

E + 2aE + a^E = a^Tr . (10) 

As a result, the model reduces to the set of equations (|6I9I10|) . For a small enough a value, the dynamics converges 
to a fixed point. Since all neurons sequentially cross the firing threshold, the translational invariance of a fixed point 
automatically implies a uniform distribution of the phase differences, that is the distinguishing feature of splay states. 
A tedious but strightforward linear stability analysis proves that the splay state destabilizes at the critical value 
already identified in ^j, where a Hopf bifurcation signals theonset of a collective periodic behaviour, i.e. of pS. 

In Fig. n we present the outcome of numerical simulations carried on in the pS regime. First of all, notice that the 
results are basically indistinguishable from those of the original set of equations, as they should. It is then instructive 
to look at T(t), since it provides a bridge with the microscopic description. In fact, if we denote with t„[i] the time 
when the ith neuron emits a spike (i.e., it crosses the threshold cf)), T{tn[i]) represents the last inter-spike interval 
(ISI) of the very same neuron and all past spike times are obtained by simply iterating the recursive formula 

i„_i[i] =t„H-T(i„[i]). (11) 

Geometrically, the implemetation of this formula corresponds to the following procedure (see Fig. ^p) : given the point 
{tn,T(tn)), move down along a straight line of slope 1 until the t-axis is reached at the time and then vertically 
up until the T-curve is hit. Considering that T{t) is a positive definite periodic function, Eq. (|ll|l describes nothing 
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FIG. 1: Comparison between the initial model for A'^ — 100 oscillators (crosses) and the set of equations 16191101 (solid line) 
for a = 0.3, g — 0.4, A = 1 and a = 9. The representation in the phase space (E,E) is plotted in panel (a), while the time 
evolution of T and n is drawn in (b) and (c). The dots in panel (a) correspond to simulations performed in the presence of 
additive uniformly distributed noise of width g = 0.025 (the wider band corresponds to A'^ = 200, while the other to A'' = 4000). 
Finally, a geometric procedure for reconstructing the single neuron firing times is plotted in (b). 



but the evolution equation of a periodically forced oscillator, so that one expects the onset of locking phenomena. 
However, it appears that the self-generated T dynamics does not exhibit any such locking, even when the ratio of the 
two frequencies is rational. This remarkable feature remains unexplained. 

Another peculiarity of pS is that the single-neuron ISI differs from (it is always smaller than) the period of the 
macroscopic dynamics. Interestingly enough, this is the distinguishing feature of slowly oscillating periodic solutions 
already found in equations with state-dependent delay (21.j . Qualitatively, the phenomenon can be clarified by noticing 
that neurons tend to bunch together, but, at the same time, those neurons lying in the front of the cluster escape 
away, while those reaching the back get stuck. This justifies the name "partial synchronization" attributed to this 
phenomenon. 



IV. PERTURBATIVE ANALYSIS 



Upon exploring the parameter plane (a. A), we have discovered that pS can be observed for arbitrarily small A- 
values. As a result, one can hope to shed further light on this phenomenon, by developing a perturbative analysis. 
This is not easy, because the limit case A = is degenerate, as Eq. ^ is satisfied by any periodic function n{t). In 
fact, in this limit, there are no effective interactions among the oscillators and their distribution does not change in 
time, whatever the initial choice is. The only constraints are given by Eq. H1U|) . which allows determining E{t) and 
Eq. lO, which fixes the period P = (1 — g)/a. Thus, exactly for A = 0, pS cannot arise, since the periodicity of the 
single oscillators necessarily coincides with the collective periodicity. However, as soon as A > 0, T can differ from 
the period P, T{t) = P + \T{t). Accordingly, for a generic function u{t) periodic of period P, one can write 

u{t - T) = u{t - At) = u{t) - \Tu{t) + 0{\^). (12) 

By introducing this expansion into the model equations H9I10II . one obtains, up to first order in A, 

1 - qrE 

T = -P+ 

a + gE 

C 

E + 2aE + a^E ^ a^ — . (13) 
r 

The same procedure, when applied to Eq. ©, implies that ttt = C, where C is a constant of motion, so that the 
model reduces to a simple set of ordinary differential equations in a three dimensional space (P, E and r being the 
relevant variables). The a priori unknown parameter C must be determined self-consistently by imposing that the 
integral oi -k = C /t over the period P is equal to 1. By numerically integrating the above set of equations, we find a 
self-consistent solution above a suitable critical a-value. The corresponding collective behaviour is plotted in Fig. |3 
for different a- values. 
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FIG. 2: The amplitude of E and r decreases as one approaches etc ~ 4.055(8) (the critical point for A — 0.001 and g = 0.4). 
The transition is thus a supercritical Hopf bifurcation. 



In order to test the robustness of the pS regime, we have performed some further studies both to verify whether the 
discontinuity in the force field is a necessary condition and to explore the stability against the addition of external 
noise. The simplest way to introduce a continuous force field is by adding a branch F'{(j)) = a + X'cf) + E (with 
(j) < = A/(A' — A)) to Eq. © which is now restricted to > (pQ. It is in principle possible to follow the same 
approach that has led to Eq. 0, but here we limit ourselves to mentioning the result of our numerical simulations, 
namely that for A' >~ 14A, the continuous model keeps exhibiting a transition to pS. Finally, we have numerically 
analysed the original model, under the further action of a zero-average uniformly distributed noise of width a — 0.025 
for N = 200 and 4000. The results, plotted in Fig. show that the collective behaviour, though depressed by the 
noise action, persists. In fact, the transversal fuzziness appears to decrease with the number N of oscillators as . 

In the second section we have obtained a closed set of equations to treat identical pulse-coupled rotators. The main 
obstacle towards specific applications is the need to express the finite-time Lyapunov exponent in Eq. (0J in terms 
of the other relevant variables. This seems to be only a technical problem, that we have indeed recently solved for 
generic sinusoidal fields |2^. On the other hand, it is less clear how to incorporate the effect of noise and disorder. 
In this respect, a close comparison with the dynamical mean field approach should be very helpful, since noise can be 
easily handled by the latter one, while our method is more open to the treatment of generic force fields. 

We thank an unknown referee for having drawn to our attention a similar bifurcation found in systems with variable 
time delay. 
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